R

R es un entorno y un lenguaje de programación estadística.

Es una de las herramientas más usadas en el análisis de datos. Es de código libre, su uso es gratuito y dispone de una comunidad muy activa que contribuye al desarrollo de herramientas y funciones adicionales (más de 10.000 paquetes) y a la resolución de preguntas a través de la red.

https://www.r-project.org/about.html https://www.rdocumentation.org/ https://stackoverflow.com/questions/tagged/r

Instalación de R y RStudio

• R

  Se puede instalar desde https://cran.r-project.org/.

• RStudio

  Es una interfaz gráfica (o IDE) para la realización de cálculos y la programación en R.

  Se instala desde https://www.rstudio.com/, versión RStudio Desktop (Open Source License).

R como herramienta de cálculo

R puede usarse desde la consola introduciendo directamente las instrucciones correspondientes.

2 + 3 

## [1] 5

8 ^ 10

## [1] 1073741824

log(100) # Logaritmo neperiano

## [1] 4.60517

# indica que lo que sigue es un comentario

Uso de scripts

Un script en R es un conjunto de instrucciones escritas en un archivo de texto, de forma que estás puedan almacenarse y ejecutarse con posterioridad. Se suelen almacenar con la extensión .R.

Un editor muy recomendado para la creación, ejecución y depuración de estos scripts es RStudio.

Acceso a la ayuda

• Buscar ayuda sobre una función

? "mean"

• Buscar un texto en la ayuda

?? "anova"

En RStudio, se puede seleccionar y pulsar F1 para buscar en la ayuda

Datos básicos

• numeric: número decimal de 15 dígitos significativos
• integer: entero de 32 bits (hasta 2·10⁹ aprox.)
• logical: TRUE o FALSE
• character: cadena de caracteres de longitud indeterminada

Datos básicos – numeric

a <- 2
a

## [1] 2

class(a)

## [1] "numeric"

b <- 13.6788956789
b

## [1] 13.6789

class(b)

## [1] "numeric"

print(b, digits = 10)

## [1] 13.67889568

Operadores más comunes para datos de tipo numeric

a+b

## [1] 15.6789

a-b

## [1] -11.6789

a*b

## [1] 27.35779

a^b

## [1] 13114.69

a/b

## [1] 0.1462106

Datos básicos – integer

Permite almacenar números enteros. Es el tipo de datos usado para contadores e índices.

n <- as.integer(340000)
class(n)

## [1] "integer"

n <- 2L
class(n)

## [1] "integer"

Datos básicos – character

Permite almacenar cadenas de texto de longitud indeterminada.

Los literales se indican entre comillas dobles (o simples).

a <- "aaa"
b <- "bbb"

paste(a, b, "hola", sep = ", ") 

## [1] "aaa, bbb, hola"

Operadores más comunes para datos de tipo character

Datos básicos – logical

Almacena valores que son verdadero (TRUE) o falso (FALSE).

Es el resultado de una comparación.

3 == 2

## [1] FALSE

b <- 3 != 2 
b

## [1] TRUE

Operadores más comunes para datos de tipo logical

Por ejemplo…

a <- TRUE
b <- F
a & b # Operador AND

## [1] FALSE

a | b # Operador OR

## [1] TRUE

!b  # Operador NOT

## [1] TRUE

Conversión a datos básicos

Las funciones de conversión de tipos en R se indican como as. seguido del tipo de dato de destino.

as.character(2)

## [1] "2"

as.numeric(TRUE)

## [1] 1

Para saber si una información o variable es de un tipo determinado puede usarse is. seguido del tipo de dato a comprobar.

Constantes

R incorpora algunas constantes preestablecidas como

pi

## [1] 3.141593

Inf
NaN
NA
NULL

También existen funciones para evaluar si un dato corresponde a una de estas constantes.

is.na(3)

## [1] FALSE

is.null(NULL)

## [1] TRUE

is.infinite(-Inf)

## [1] TRUE

Datos compuestos

• Principales
  • Vector
  • Factor
  • Data frame (Tabla de datos)
• Otros tipos comunes
  • Factor ordenado
  • Lista
  • Matriz

Datos compuestos – vector

a <- c(2, 3, 4)
str(a)

##  num [1:3] 2 3 4

a[2]

## [1] 3

Las operaciones se aplican a vectores y devuelven vectores

a + a 

## [1] 4 6 8

a > 2.5

## [1] FALSE  TRUE  TRUE

Atención al reciclaje de datos…

a <- c(2, 3, 4)
b <- c(10, 20)
a * b

## Warning in a * b: longer object length is not a multiple of shorter object
## length

## [1] 20 60 40

Las funciones en R suelen ser vectoriales

• algunas devuelven vectores

a <- c(2, 3, 4)
abs(sin(a))

## [1] 0.9092974 0.1411200 0.7568025

exp(a)

## [1]  7.389056 20.085537 54.598150

• otras devuelven valores agregados

length(a)

## [1] 3

sum(a)

## [1] 9

mean(a)

## [1] 3

También sd, max, min…

Puede determinarse si un valor está presente en un vector con el operador %in%.

a <- c(2, 3, 4)
3 %in% a

## [1] TRUE

c(2,5,3,4) %in% a

## [1]  TRUE FALSE  TRUE  TRUE

Existen estructuras especificas para crear vectores secuencia.

1:10

##  [1]  1  2  3  4  5  6  7  8  9 10

seq(1, 10, by = 2)

## [1] 1 3 5 7 9

seq(1, 3, length.out = 5) 

## [1] 1.0 1.5 2.0 2.5 3.0

Los vectores se ordenan usando las funciones sort y order.

a <- c(8, 2, 5, 3)
sort(a)

## [1] 2 3 5 8

a[order(a)]

## [1] 2 3 5 8

a[order(a,decreasing = T)]

## [1] 8 5 3 2

Datos compuestos – factor

Un factor es un vector de cadenas indexado. Normalmente se crea a partir del vector de cadenas de texto.

a <- c("hola", "adeu","hola", "adeu", "adeu", "bye")
b <- as.factor(a)
as.character(b)

## [1] "hola" "adeu" "hola" "adeu" "adeu" "bye"

as.numeric(b)

## [1] 3 1 3 1 1 2

str(b)

##  Factor w/ 3 levels "adeu","bye","hola": 3 1 3 1 1 2

La función factor permite la codificación o recodificación manual de un vector.

a <- factor(c(3, 1, 3, 1, 1, 2), labels = c("adeu", "bye", "hola"))
a

## [1] hola adeu hola adeu adeu bye 
## Levels: adeu bye hola

levels(a)

## [1] "adeu" "bye"  "hola"

Datos compuestos – data frame

dfA <- data.frame(int = 1:10,
                  let = sample(letters, 10, replace = TRUE), 
                  ran = rnorm(10))
dfA

##    int let         ran
## 1    1   t -1.72203721
## 2    2   s -1.92532196
## 3    3   r -0.27530247
## 4    4   z  0.46023366
## 5    5   a  0.56863506
## 6    6   e  1.28603975
## 7    7   p -1.56824221
## 8    8   w  0.01242195
## 9    9   g  1.05118640
## 10  10   b  1.28818348

dim(dfA) # Dimensión

## [1] 10  3

nrow(dfA) # Número de columnas

## [1] 10

ncol(dfA) # Número de filas

## [1] 3

str(dfA)

## 'data.frame':    10 obs. of  3 variables:
##  $ int: int  1 2 3 4 5 6 7 8 9 10
##  $ let: Factor w/ 10 levels "a","b","e","g",..: 8 7 6 10 1 3 5 9 4 2
##  $ ran: num  -1.722 -1.925 -0.275 0.46 0.569 ...

head(dfA, 3) # Primeros datos, 6 por defecto

##   int let        ran
## 1   1   t -1.7220372
## 2   2   s -1.9253220
## 3   3   r -0.2753025

tail(dfA, 2) # Últimos valores

##    int let      ran
## 9    9   g 1.051186
## 10  10   b 1.288183

Para acceder a los datos almacenados en el data frame se usan índices. Las variables de data frame también pueden extraerse usando su nombre.

dfA[2,3]

## [1] -1.925322

dfA[,1]

##  [1]  1  2  3  4  5  6  7  8  9 10

dfA$let

##  [1] t s r z a e p w g b
## Levels: a b e g p r s t w z

Veremos más formas de acceder y usar los datos de una tabla de datos más adelante.

Más datos compuestos – lista

a <- list(2, "2", FALSE)
b <- list(3, "hola", c(2, 3, 4))

a

## [[1]]
## [1] 2
## 
## [[2]]
## [1] "2"
## 
## [[3]]
## [1] FALSE

b

## [[1]]
## [1] 3
## 
## [[2]]
## [1] "hola"
## 
## [[3]]
## [1] 2 3 4

Más datos compuestos – factor ordenado

notes <- c("Aprovat", "Insuficient", "Notable", "Insuficient",
           "Notable", "Excel·lent", "Aprovat")
notes <- factor(notes,
        levels = c("Insuficient", "Aprovat",
                   "Notable", "Excel·lent"),
        ordered = TRUE)
str(notes)

##  Ord.factor w/ 4 levels "Insuficient"<..: 2 1 3 1 3 4 2

levels(notes)

## [1] "Insuficient" "Aprovat"     "Notable"     "Excel·lent"

Más datos compuestos – matriz

a <- matrix(c(2, 4, -3, 5), ncol = 2)
a

##      [,1] [,2]
## [1,]    2   -3
## [2,]    4    5

a[2,2]

## [1] 5

a * a # Producto posición por posición

##      [,1] [,2]
## [1,]    4    9
## [2,]   16   25

a %*% a # Producto matricial

##      [,1] [,2]
## [1,]   -8  -21
## [2,]   28   13

t(a) # Transposición

##      [,1] [,2]
## [1,]    2    4
## [2,]   -3    5

Instalación y carga de paquetes en R

R tiene muchos paquetes para resolver problemas específicos. Para usar un paquete adicional este debe instalarse y cargarse en memoria.

installed.packages()[,1] # Lista los paquetes instalados
(.packages()) # Lista los paquetes en memoria

Dos recursos importantes para buscar y identificar paquetes relevantes son https://www.rdocumentation.org/ https://cran.rstudio.com/web/views/

Para instalar un paquete, por ejemplo “nycflights13”

install.packages("nycflights13")

Para cargar un paquete en memoria, se usa

library("nycflights13")

En RStudio la gestión de paquetes también puede hacerse desde la interfaz del programa

En un script y para garantizar la disponibilidad de un paquete

if(!require("nycflights13")) {
  install.packages("nycflights13")
  library("nycflights13")
}

Para acceder a la documentación de un paquete

help(package="nycflights13")

Algunos paquetes tienen información adicional a la que se puede acceder con las funciones vignette() y demo().

¿Para qué usamos los gráficos?

La principales funciones de los gráficos en el análisis de datos son

• explorar los datos y facilitar su comprensión,
• permitir el descubrimiento de patrones no evidentes, y
• comunicar efectivamente los resultados de los análisis.

Gráficos en R

Veremos dos paradigmas distintos para la creación de gráficos en R.

• Gráficos creados con el paquete base.
• Gráficos creados en base a la gramática de gráficos (GoG), usando el paquete ggplot2

Gráficos en base R

Los gráficos en base R se construyen a partir de datos en vectores y usando funciones específicas en función del gráfico a realizar.

Para ejemplificar algunas de estas funciones, usaremos el conjunto de datos mtcars, que forma de los paquetes básicos de R.

str(mtcars)

## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

head(mtcars, 10)

##                    mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360        14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 240D         24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230          22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280          19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4

Veremos como crear

• gráficos de dispersión,
• histogramas,
• diagramas de barras, y
• diagramas de caja (boxplot)

Gráficos en base R – gráfico de dispersión

plot(mtcars$disp,mtcars$hp)

Gráficos en base R – histograma

hist(mtcars$disp)

Añadiendo una curva de densidad…

hist(mtcars$disp,breaks=20,freq=FALSE)
lines(density(mtcars$disp))

Gráficos en base R – diagrama de barras

barplot(table(as.factor(mtcars$cyl)))

Gráficos en base R – diagrama de caja

boxplot(mtcars$disp)
points(mean(mtcars$disp))

Gramática de gráficos (GoG)

La gramática de gráficos es una aproximación teórica al estudio de los componentes de un gráfico. De acuerdo con este análisis, un gráfico se puede construir mediante la especificación de un conjunto de capas y componentes que definen los datos, la asociación de los mismos a aspectos del gráfico, la especificación de la relación entre los valores de las variables de los datos con las del gráfico, la estructura geométrica del gráfico…

Wilkinson, L. (2006). The grammar of graphics. Springer Science & Business Media.

ggplot2

ggplot2 es un paquete de R que implementa de la gramática de gráficos.

Wickham, H. (2010). A layered grammar of graphics. Journal of Computational and Graphical Statistics, 19(1), 3-28.

Referencias:

• http://ggplot2.tidyverse.org/reference/
• https://github.com/rstudio/cheatsheets/raw/master/data-visualization-2.1.pdf
• http://www.ling.upenn.edu/~joseff/avml2012/
• http://r-statistics.co/Complete-Ggplot2-Tutorial-Part1-With-R-Code.html

ggplot2 forma parte del paquete tidyverse (aunque también puede instalarse y cargarse autónomamente).

if(!require("tidyverse")) {
  install.packages("tidyverse")
  library("tidyverse")
}

ggplot2 – capas y elementos del gráfico

En ggplot2 cada elemento gráfico que representa un conjunto de datos constituye una capa. Una o más capas constituyen un gráfico.

Cada capa queda definida mediante la especificación de sus elementos. Los principales son

• datos,
• mapeado estético, y
• geometrías

En ggplot2, los gráficos constituyen un objeto de R y se construyen de forma aditiva.

Por ejemplo,

grafico <- ggplot(data = anscombe,
        mapping = aes(x = x1, y = y1))  # Datos y mapeado estético 
grafico <- grafico + geom_point()       # Geometría

grafico

ggplot2 – datos

En ggplot2, el elemento data (datos) se introduce como primer argumento de la función ggplot. Debe corresponder a una tabla de datos o un tipo de datos convertible a tabla de datos.

grafico <- ggplot(data = anscombe,

ggplot2 – mapeado estético

El mapping (mapeado estético) corresponde al establecimiento de relaciones entre variables de los datos y variables del gráfico. Es el segundo argumento de la función ggploty debe crearse con al función de apoyo aes.

        mapping = aes(x = x1, y = y1))

Para variables cuantitativas, los mapeados más comunes corresponden a

• posiciones: x, y…
• tamaño: size
• color: color, fill

Para variable cualitativas, los mapeados más frecuentes son

• posiciones: x, y…
• color: color, fill
• forma: shape

ggplot2 – geometrías

Las geometrías (geom_) indican la forma que debe tener el gráfico, es decir, cómo se articulan las variables del gráfico. Se añaden al gráfico sumándose al objeto creado por ggplot.

grafico <- grafico + geom_point()

Son geometrías de uso común

• geom_point
• geom_line, geom_vline, geom_hline
• geom_bar
• geom_histogram
• geom_boxplot

Un resumen de las geometrías y su relación con las variables del gráfico que reconoce cada una de ellas figura en https://github.com/rstudio/cheatsheets/raw/master/data-visualization-2.1.pdf.

Gráficos en ggplot2

Veremos cómo crear en ggplot2 los gráficos más habituales, añadiendo algunas consideraciones para aquellos casos donde los gráficos realizados con base tienen prestaciones insuficientes.

• gráficos de dispersión,
• histogramas,
• diagramas de barras, y
• diagramas de caja (boxplot).

De forma general, los gráficos en ggplot2 se construyen a partir de tablas de datos (data frames), de los cuales se seleccionan las variables a representar.

Usaremos 1000 datos del conjunto de datos diamonds para crear los distintos ejemplos.

diaM <- diamonds[sample(1:nrow(diamonds),1000),]
str(diaM)

## Classes 'tbl_df', 'tbl' and 'data.frame':    1000 obs. of  10 variables:
##  $ carat  : num  0.54 0.4 2.21 2.4 0.55 0.33 1.5 2.63 0.52 0.58 ...
##  $ cut    : Ord.factor w/ 5 levels "Fair"<"Good"<..: 4 3 5 5 5 5 3 5 5 5 ...
##  $ color  : Ord.factor w/ 7 levels "D"<"E"<"F"<"G"<..: 6 6 6 6 1 4 2 7 4 3 ...
##  $ clarity: Ord.factor w/ 8 levels "I1"<"SI2"<"SI1"<..: 6 4 3 2 3 6 3 2 7 4 ...
##  $ depth  : num  60.9 61.9 61.6 63.8 60.9 61.3 58.6 62.1 61.6 62.1 ...
##  $ table  : num  59 55 57 58 57 56 62 57 56 55 ...
##  $ price  : int  1944 1031 18062 17039 1840 730 10626 16914 1881 1790 ...
##  $ x      : num  5.28 4.73 8.38 8.41 5.29 4.46 7.43 8.83 5.16 5.37 ...
##  $ y      : num  5.26 4.77 8.43 8.37 5.32 4.48 7.46 8.88 5.2 5.33 ...
##  $ z      : num  3.21 2.94 5.18 5.35 3.23 2.74 4.36 5.5 3.19 3.32 ...

Gráficos en ggplot2 – gráfico de dispersión

ggplot(diaM, aes(x=carat,y=price)) + geom_point()

Añadiendo una tercera variable (cut) y modificando algunos aspectos de formato…

ggplot(diaM, aes(x=carat,y=price,color=cut)) + 
  geom_point(alpha=.8,shape=21,size=3)

Añadiendo líneas de tendencia…

ggplot(diaM, aes(x=carat,y=price,color=cut)) + 
  geom_point(alpha=.8,shape=21,size=3) +
  geom_smooth(method="lm",se=FALSE)

Gráficos en ggplot2 – histograma

ggplot(diaM, aes(x=price)) + geom_histogram(binwidth=1000)

Y en función del corte…

ggplot(diaM, aes(x=price,fill=cut)) +
  geom_histogram(position='dodge',binwidth=1000)

En frecuencias relativas (por grupo)…

ggplot(diaM, aes(x=price,y=..density..,fill=cut)) +
  geom_histogram(position='dodge',binwidth=1000)

Quizás funcione mejor un gráfico de densidades…

ggplot(diaM, aes(x=price,fill=cut)) +
  geom_density(alpha=.3)

Gráficos en ggplot2 – diagrama de barras

ggplot(diaM, aes(x=clarity)) + geom_bar()

En función de la claridad…

ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar()

Para comparar entre frecuencias absolutas, funcionan mejor las barras separadas.

ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar(position="dodge")

Para comparar entre frecuencias relativas acumuladas, son mejores las barras apiladas en frecuencia relativa (para cada clase).

ggplot(diaM, aes(x=clarity, fill=cut)) + geom_bar(position="fill")

Los diagramas de barras también pueden crearse a partir de tablas de datos agrupados. En este caso, debe indicarse qué variable es la y e incluir stat="identity" en el geom_bar.

Gráficos en ggplot2 – diagrama de caja

ggplot(diaM, aes(x=1, y=price)) + geom_boxplot()

NOTA: Para el geom_boxplot se requieren dos variables. Si no hay variable independiente se puede incluir x=1.

Y en función del corte…

ggplot(diaM, aes(x=cut, y=price)) + geom_boxplot()

El gráfico se puede mejorar mostrando todos los puntos, con una posición aleatorizada.

ggplot(diaM, aes(x=cut, y=price)) + 
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(shape = 21, alpha=.5,height=0,width=.2)

O incluyendo un violin plot y un punto para la media…

medias <- diaM %>% group_by(cut) %>%
  summarise(price=mean(price))

ggplot(diaM, aes(x=cut, y=price)) + 
  geom_violin() +
  geom_boxplot(outlier.shape = NA, width = 0.1) +
  geom_point(data=medias,shape=3)

ggplot2 – elementos adicionales

Además de los elementos ya presentados (datos, mapeado estético y geometrías), otros elementos de ggplot2 permiten controlar aspectos adicionales del gráfico, por ejemplo

• scale_...: controlan los aspectos relativos a la presentación de las variables del gráfico,
• coord_...: establecen el sistema de coordenadas usado en la geometría,
• labs: establece los títulos del gráfico,
• theme, theme_...: controlan la part del gráfico que no corresponde a datos (non-data ink), y
• facet_: permite la creació de secuencia de gráficos en función de una o dos variables

Un ejemplo para terminar…

ggplot(diaM, aes(x=carat, y = price, shape = cut, col = clarity)) +
  geom_point(alpha=.6) +
  scale_x_continuous(breaks=1:3) +
  scale_y_continuous(trans="log10") +
  scale_color_brewer(palette="Spectral")+
  facet_grid(cut ~ clarity) + 
  theme_bw() + 
  theme(legend.position = "none",text = element_text(size=10))

Funciones de la estadística

La estadística, como parte de la matemática que se ocupa de la recolección, análisis e interpretación de datos, tiene dos funciones principales

1. la descripción de conjuntos de datos (estadística descriptiva), y
2. la extracción de conclusiones a partir de estos datos (inferencia)

En este repaso de la estadística básica, usaremos (de nuevo) el conjunto de datos mtcars.

str(mtcars)

## 'data.frame':    32 obs. of  11 variables:
##  $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
##  $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
##  $ disp: num  160 160 108 258 360 ...
##  $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
##  $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
##  $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
##  $ qsec: num  16.5 17 18.6 19.4 17 ...
##  $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
##  $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
##  $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
##  $ carb: num  4 4 1 1 2 1 4 2 2 4 ...

head(mtcars, 10)

##                    mpg cyl  disp  hp drat    wt  qsec vs am gear carb
## Mazda RX4         21.0   6 160.0 110 3.90 2.620 16.46  0  1    4    4
## Mazda RX4 Wag     21.0   6 160.0 110 3.90 2.875 17.02  0  1    4    4
## Datsun 710        22.8   4 108.0  93 3.85 2.320 18.61  1  1    4    1
## Hornet 4 Drive    21.4   6 258.0 110 3.08 3.215 19.44  1  0    3    1
## Hornet Sportabout 18.7   8 360.0 175 3.15 3.440 17.02  0  0    3    2
## Valiant           18.1   6 225.0 105 2.76 3.460 20.22  1  0    3    1
## Duster 360        14.3   8 360.0 245 3.21 3.570 15.84  0  0    3    4
## Merc 240D         24.4   4 146.7  62 3.69 3.190 20.00  1  0    4    2
## Merc 230          22.8   4 140.8  95 3.92 3.150 22.90  1  0    4    2
## Merc 280          19.2   6 167.6 123 3.92 3.440 18.30  1  0    4    4

Descriptiva de una variable

A menudo, el primer interés que tendremos es el de describir un conjunto de datos correspondientes a una sola variable. Esta es la situación con la que empezaremos.

¿Cómo se describe un conjunto de datos univariante?

Normalmente se consideran los siguientes aspectos

• frecuencias y cuantiles,
• tendencia central (también llamado localización o posición),
• dispersión, y
• análisis de puntos extremos

Descriptiva de una variable – frecuencias y cuantiles

Tamaño de la muestra o número de datos

length(mtcars$disp)

## [1] 32

n <- sum(!is.na(mtcars$disp))
n

## [1] 32

Frecuencias absolutas y relativas

Para una variable cuantitativa (numeric)…

frec_abs <- table(cut(mtcars$disp,
      breaks=c(0,100,200,300,400,500,Inf)))
frec_abs

## 
##   (0,100] (100,200] (200,300] (300,400] (400,500] (500,Inf] 
##         5        11         5         8         3         0

frec_rel <- frec_abs / n
frec_rel

## 
##   (0,100] (100,200] (200,300] (300,400] (400,500] (500,Inf] 
##   0.15625   0.34375   0.15625   0.25000   0.09375   0.00000

Cuantiles

quantile(mtcars$disp,.1)

##   10% 
## 80.61

quantile(mtcars$disp,1:4*.2)

##    20%    40%    60%    80% 
## 120.14 160.00 275.80 350.80

Descriptiva de una variable – tendencia central

mean(mtcars$disp)

## [1] 230.7219

median(mtcars$disp)

## [1] 196.3

Descriptiva de una variable – dispersión

var(mtcars$disp) # el denominador es n-1

## [1] 15360.8

sd(mtcars$disp)

## [1] 123.9387

range(mtcars$disp)

## [1]  71.1 472.0

diff(range(mtcars$disp))

## [1] 400.9

IQR(mtcars$disp)

## [1] 205.175

mad(mtcars$disp) # mean absolute deviation

## [1] 140.4764

Descriptiva de una variable – puntos extremos

lmin <- quantile(mtcars$disp,.25) - 1.5*IQR(mtcars$disp) 
lmax <- quantile(mtcars$disp,.75) + 1.5*IQR(mtcars$disp) 
mtcars$disp[mtcars$disp > lmax | mtcars$disp < lmin]

## numeric(0)

lmin <- quantile(mtcars$qsec,.25) - 1.5*IQR(mtcars$qsec) 
lmax <- quantile(mtcars$qsec,.75) + 1.5*IQR(mtcars$qsec) 
mtcars$qsec[mtcars$qsec > lmax | mtcars$qsec < lmin]

## [1] 22.9

Descriptiva de una variable – gráficos

Gráficamente, son útiles para describir un conjunto de datos

• para variables cuantitativas: el diagrama de caja y el histograma

boxplot(mtcars$qsec)
hist(mtcars$qsec)

• para variables cualitativas: el diagrama de barras

barplot(table(factor(mtcars$cyl)))

Descriptiva de dos variables

Para el estudio y descripción de la relación entre dos variables, las técnicas más habituales incluyen

• las tablas de contingencia,
• los coeficientes de correlación, y muy especialmente el producto-momento de Pearson, y
• los gráficos de dispersión

Descriptiva de dos variables – tabla de contingencia

Es el concepto equivalente a las tablas de frecuencias en el caso univariante.

Para variables cualitativas…

table(factor(mtcars$am),factor(mtcars$cyl))

##    
##      4  6  8
##   0  3  4 12
##   1  8  3  2

Para variables cuantitativas…

table(cut(mtcars$disp,breaks = seq(0,500,by=100)),
      cut(mtcars$qsec,breaks = seq(13,25,by=3)))

##            
##             (13,16] (16,19] (19,22] (22,25]
##   (0,100]         0       3       2       0
##   (100,200]       1       7       2       1
##   (200,300]       0       3       2       0
##   (300,400]       4       4       0       0
##   (400,500]       0       3       0       0

Descriptiva de dos variables – coeficiente de correlación

Producto-momento de Pearson

cor(mtcars$disp,mtcars$qsec)

## [1] -0.4336979

Coeficiente de correlación de Spearman

cor(mtcars$disp,mtcars$qsec, method = "spearman")

## [1] -0.4597818

Descriptiva de dos variables – gráfico de dispersión

plot(mtcars$disp,mtcars$qsec)

Distribuciones de probabilidad

La distribución de probabilidad de un número aleatorio corresponde a la abstracción teórica de la densidad de los datos de una conjunto. Representa la densidad que se obtendría cunado se generasen infinitos números aleatorios de acuerdos a un mismo procedimiento o experimento.

Un gran número de distribuciones experimentales tienen modelos teóricos con los que se relacionan. Los tres más comunes son la distribución normal (para datos de variables cuantitativas continuas), la distribución uniforme (para cualquier tipo de datos) y la distribución binomial (para el número de veces que se produce un suceso de una determinada probabilidad).

En R, todas las distribuiciones comparten el mismo sistema de funciones

• r<dist>: genera aleatorios deacuerdo con una distribución de probabilidad
• d<dist>: devuelve la densidad de probablidad para un valor de variable
• q<dist>: devuelve el valor de la variable tal que la probablidad
• p<dist>: devuelve la probablidad que un valor sea inferior o igual a x

Por ejemplo, para la distribución normal estándar…

rnorm(10)

##  [1] -1.4887791  0.2170026  0.1122813 -1.5607699  1.2536196 -0.4178993
##  [7] -0.9548049 -0.6038304 -2.3527068  0.8741368

dnorm(0)

## [1] 0.3989423

qnorm(.95)

## [1] 1.644854

pnorm(1.64)

## [1] 0.9494974

Distribuciones de probabilidad - normal

df <- data.frame(x = rnorm(1000, mean = 3, sd = 1))
dfT <-data.frame(x = seq(0,6,length.out=101),
      y = dnorm(seq(0,6,length.out=101),mean=3,sd=1))
ggplot(df, aes(x = x, y=..density..)) +
  geom_histogram(binwidth = .1,fill="#dddddd",col="black") +
  geom_density() +
  geom_line(data=dfT,aes(x=x,y=y),col="red")+
  theme_bw()

pnorm(5,mean = 3,sd = 1)

## [1] 0.9772499

qnorm(.98,mean = 3,sd = 1)

## [1] 5.053749

Distribuciones de probabilidad - uniforme

df <- data.frame(x = runif(1000, min = 10, max = 20))
dfT <-data.frame(x = seq(10,20,length.out=101),
      y = dunif(seq(10,20,length.out=101), min=10, max=20))
ggplot(df, aes(x = x, y=..density..)) +
  geom_histogram(binwidth = .5,fill="#dddddd",col="black") +
  geom_density() +
  geom_line(data=dfT,aes(x=x,y=y),col="red")+
  theme_bw()

punif(12, min=10, max=20)

## [1] 0.2

qunif(.90, min=10, max=20)

## [1] 19

Distribuciones de probabilidad - binomial

df <- data.frame(x = rbinom(100,5,prob=0.5))
df <- df %>% group_by(x) %>%
  summarise(y = n()/nrow(df))
dfT <- data.frame(x = 0:5,
    y = dbinom(0:5,size=5,prob=0.5))

dfJ <- rbind(cbind(teo=FALSE,df),
             cbind(teo=TRUE,dfT))

ggplot(dfJ,aes(x=x,y=y,fill=teo)) +
  geom_bar(position="dodge", stat="identity",
           col="black") +
  scale_fill_manual(values=c("#dddddd","red")) +
  theme_bw() +
  theme(legend.position = "none")

pbinom(2,5,prob=0.5)

## [1] 0.5

qbinom(.5,5,.5)

## [1] 2

Distribuciones de probabilidad – otras

R incorpora otras muchas distribuciones de probabilidad que pueden consultarse en la página correspondiente de la ayuda.

? "Distributions"

Inferencia

En términos estadísticos, la inferencia consiste en la obtención de información sobre la población –el conjunto de los valores existentes– a partir de una muestra representativa de la misma.

Si los datos que tenemos, no constituyen una muestra representativa cualquier información que queramos extrapolar a la población estará sesgada, es decir, desplazada de su valor real.

Los resultados de la inferencia se suelen mostrar mediante una de las dos formas siguientes:

• el valor de un parámetro para la población, asociado normalmente a un error o un intervalo de confianza,
• la probabilidad (p) que unos datos, o cualquiera más extremo que ellos, se produzcan siendo cierta una determinada hipótesis, denominada hipótesis nula.

Presentaremos a continuación, las técnicas más comúnmente usadas para

• analizar el ajuste a una distribución,
• estudiar la dispersión de la población, y
• estimar o comparar valores de tendencia central –localización o posición–

Inferencia – distribución

Para contrastar la hipótesis que un conjunto de datos puede representarse mediante un determinada distribución teórica, las pruebas más comunes son la prueba de bondad de ajuste de chi cuadrado (chisq.test) y la prueba de Kolmogorov-Smirnov (ks.test). En ambos casos, los argumentos son las frecuencias absolutas de cada clase para la distribución experimental y las probabilidades esperadas de acuerdo con la distribución teórica o de referencia.

Ejemplo

Se ha tirado cien veces un dado de 10 caras, y se han obtenido cada una de las caras las veces que se muestran en la tabla siguiente.

Si el dado fuese normal, todas las caras tendrían la misma probabilidad de aparecer. ¿Lo es?

obs <- c(9, 6, 7, 15, 11, 11, 9, 13, 9, 10)  
exp <- rep(.1, 10)
 
chisq.test(x = obs, p = exp)

## 
##  Chi-squared test for given probabilities
## 
## data:  obs
## X-squared = 6.4, df = 9, p-value = 0.6993

Gráficamente, y aunque sin la misma objetividad, la comparación de distribuciones también suele hacerse mediante el gráfico de cuantiles –en R, es la función qqplot.

Si ambos conjuntos de datos corresponden a la misma distribución, sus cuantiles deben quedar alineados junto a la diagonal del gráfico (especialemente en la región central del mismo).

exp <- rep(1:10, obs)
teo <- rep(1:10, 10)
qqplot(x = exp, y = quantile(teo, (1:length(exp))/length(exp)))
qqline(y = exp, distribution = function(x)  quantile(teo,x),
       probs= c(1/length(exp),1))

Para comprobar si un conjunto de datos puede proceder de una población aleatoria de distribución normal, existen pruebas estadísticas más específicas. Entre las muchas pruebas existentes, son de uso común las pruebas de Shapiro-Wilk y Anderson-Darling.

x1 <- rnorm(20, mean = 4, sd = 5)
x2 <- rbeta(20, shape1 = 5, shape2 = .5, ncp = 4)

shapiro.test(x1)

## 
##  Shapiro-Wilk normality test
## 
## data:  x1
## W = 0.96718, p-value = 0.6945

shapiro.test(x2)

## 
##  Shapiro-Wilk normality test
## 
## data:  x2
## W = 0.84082, p-value = 0.003755

if(!require("nortest")) {
  install.packages("nortest")
  library("nortest")
}

ad.test(x1)

## 
##  Anderson-Darling normality test
## 
## data:  x1
## A = 0.25847, p-value = 0.6791

ad.test(x2)

## 
##  Anderson-Darling normality test
## 
## data:  x2
## A = 1.2518, p-value = 0.002178

Gráficamente se puede usar el gráfico de cuantiles comentado anteriormente.

qqnorm(x1)
qqline(x1)

qqnorm(x2)
qqline(x2)

Inferencia – dispersión

Las inferencias más habituales respecto a la dispersión de una población constituyen los siguientes casos

• determinar el intervalo de confianza de la varianza o la desviación típica de la distribución,
• comparar la dispersión de dos conjuntos de datos,
• comparar la dispersión de tres o más conjuntos de datos –o comprobar la homocedasticidad de distintos conjuntos de datos–.

Cuando los datos estan normalmente distribuidos, el intervalo de confianza de la varianza se calcula a partir de la distribución de chi cuadrado.

x <- rnorm(50, mean = 0, sd = 2)
df <- length(x) - 1
lower <- var(x) * df / qchisq(1 - 0.05/2, df)
upper <- var(x) * df / qchisq(0.05/2, df)
c(lower = lower, var = var(x), upper = upper)

##    lower      var    upper 
## 3.621426 5.189904 8.059134

Para obtener el intervalo de confianza para la desviación estándar, cuando los datos estan normalmente distribuidos, se calcula la desviación estándar a partir de las varainzas calculadas más arriba.

c(lower = sqrt(lower), sd = sd(x), upper = sqrt(upper))

##    lower       sd    upper 
## 1.903005 2.278136 2.838861

Para la comparación de la dispersión de las poblaciones para dos conjuntos de datos, cuando los datos estan normalmente distribuidos, se lleva a cabo con un prueba de F –función var.test en R–.

x <- rnorm(50, mean = 0, sd = 2)
y <- rnorm(30, mean = 1, sd = 1)
var.test(x, y)

## 
##  F test to compare two variances
## 
## data:  x and y
## F = 3.9057, num df = 49, denom df = 29, p-value = 0.0001963
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  1.962337 7.348335
## sample estimates:
## ratio of variances 
##           3.905746

Si los datos no estan normalmente distribuidos, la comparación de dispersiones puede hacerse mediante el test de Ansari.

x <- rlnorm(50, meanlog = 2, sdlog = 1)
y <- rlnorm(30, meanlog = 2, sdlog = .2)
ansari.test(x, y)

## 
##  Ansari-Bradley test
## 
## data:  x and y
## AB = 746, p-value = 2.91e-08
## alternative hypothesis: true ratio of scales is not equal to 1

Para comparar las dispersiones de más de un conjunto de datos para los cuáles se pueda asumir normalidad, se usa la prueba de Bartlett –bartlett.test en R–.

x1 <- round(rnorm(20, mean = 1, sd = 2),1)
x2 <- round(rnorm(20, mean = 3, sd = 2),1)
x3 <- round(rnorm(20, mean = 5, sd = 2),1)
list(x1,x2,x3)
bartlett.test(list(x1,x2,x3))

## [[1]]
##  [1]  0.2 -2.0 -0.5  4.8 -1.9  1.1  0.5 -2.3  0.2  3.1 -1.7 -0.3  1.1  0.7
## [15] -0.9  0.7  3.0  2.7  1.7 -2.3
## 
## [[2]]
##  [1]  2.9  1.3  0.7  2.6  5.0  1.3  1.0  1.3  2.7  3.9 -0.3  3.0  1.7  0.6
## [15]  4.4  1.0 -3.0 -0.4 -0.3  3.4
## 
## [[3]]
##  [1] 6.3 4.8 7.3 3.2 1.3 4.0 2.3 5.2 7.6 4.5 6.5 7.6 3.1 5.4 6.0 2.6 5.0
## [18] 6.8 5.3 5.4

## 
##  Bartlett test of homogeneity of variances
## 
## data:  list(x1, x2, x3)
## Bartlett's K-squared = 0.15182, df = 2, p-value = 0.9269

Si los datos no cumplen con la hipótesis de normalidad, en lugar de la prueba de Bartlett se usa la prueba de Levene o la de Brown–Forsythe –leveneTest en el paquete car de R–.

Gráficamente y aunque de forma menos objetiva, la comparación de dispersiones puede hacerse usando gráficos de caja, bien directamente, bien centrando los datos.

x1 <- round(rnorm(20, mean = 1, sd = 2),1)
x2 <- round(rnorm(20, mean = 3, sd = 2),1)
x3 <- round(rnorm(20, mean = 5, sd = 2),1)

xs <- data.frame(
  group = factor(c(rep(1, length(x1)),rep(2, length(x2)),
            rep(3, length(x3)))),
  y = c(x1,x2,x3),
  centered = c(x1 - median(x1),x2 - median(x2),x3 - median(x3)))

ggplot(xs, aes(x=group,y=y)) + geom_boxplot() +
  coord_flip() + theme_bw()
ggplot(xs, aes(x=group,y=centered)) + geom_boxplot() +
  coord_flip() + theme_bw()

Inferencia – localización

Las principales inferencias respecto a la localización –o tendencia central– de los datos corresponden a

• establecer un intervalo de confianza para la tendencia central de la distribución,
• comparar la tendencia central de un conjunto de datos con un valor preestablecido,
• comparar la tendencia central de dos conjuntos de datos,
• comparar la tendencia central de tres o más conjuntos de datos.

Para establecer el intervalo de confianza de la media, cuando los datos están normalmente distribuidos, se usa la distribución t. El intervalo puede calcularse a partir de la función qt o visualizarlo como resultado de la función t.test.

x <- rnorm(10, mean = 2, sd = .5)
x

##  [1] 2.4881314 1.7967923 1.9809723 2.6532510 1.9285462 1.4351935 0.6875142
##  [8] 2.4288680 2.6855859 1.6888404

t.test(x)

## 
##  One Sample t-test
## 
## data:  x
## t = 10.028, df = 9, p-value = 3.495e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  1.531327 2.423412
## sample estimates:
## mean of x 
##   1.97737

Si los datos no están normalmente distribuidos, la función wilcox.test ofrece entre sus resultados el intervalo de confianza para la mediana.

x <- rnorm(10, mean = 3, sd = .5) ^ 3
x

##  [1] 19.429765 39.367781  8.161071 28.128441 48.464746 13.157887 16.326398
##  [8] 27.739158 33.964697 14.385759

wilcox.test(x, conf.int = TRUE)

## 
##  Wilcoxon signed rank test
## 
## data:  x
## V = 55, p-value = 0.001953
## alternative hypothesis: true location is not equal to 0
## 95 percent confidence interval:
##  14.74214 33.94726
## sample estimates:
## (pseudo)median 
##        23.7791

Las dos funciones usadas para obtener los intervalos de confianza permiten comparar la tendencia central de un conjunto de datos con un valor preestablecido. Como se ha indicado antes, la prueba de t require normalidad de los datos; la prueba de Wilcoxon, no.

x <- rnorm(10, mean = 2, sd = .5)
x

##  [1] 1.678513 1.282417 2.391797 1.868252 2.358359 2.328773 2.560042
##  [8] 2.180218 1.567101 1.663925

t.test(x, mu = 1.5)

## 
##  One Sample t-test
## 
## data:  x
## t = 3.5813, df = 9, p-value = 0.005919
## alternative hypothesis: true mean is not equal to 1.5
## 95 percent confidence interval:
##  1.679729 2.296150
## sample estimates:
## mean of x 
##   1.98794

wilcox.test(x ^ 3, mu = 8)

## 
##  Wilcoxon signed rank test
## 
## data:  x^3
## V = 33, p-value = 0.625
## alternative hypothesis: true location is not equal to 8

Para comparar la localización de dos conjuntos de datos, deben tenerse en cuenta los siguientes aspectos

• ¿los datos de ambas muestras están asociados?
• ¿los datos estan normalmente distribuidos?
• ¿tienen ambos conjuntos de datos la misma dispersión?

Si los datos están asociados, es decir que existen pares de datos en ambas muestras tales que comparten fuentes de variación, entonces

• si ambas muestras están normalmente distribuidas, o la diferencia entre conjuntos de datos está normalmente distribuido, se usa la prueba de t sobre la diferencia de los valores entre ambos conjuntos de datos. En R, puede calcularse la diferencia y luego aplicar la función t.test o bien usar directamente la función t.test con la opción paired = TRUE,
• si los datos no están normalmente distribuidos se usa la prueba de Wilcoxon –wilcox.test sobre la resta o con la opción paired = TRUE en R–.

x <- rnorm(10, mean = 3, sd = 1)
y <- rnorm(10, mean = 3.5, sd = 1)
list(x = x,y = y)

## $x
##  [1] 2.270545 2.585299 2.303677 3.348109 2.966622 3.697185 4.330774
##  [8] 3.153001 4.234848 3.071585
## 
## $y
##  [1] 3.437685 2.001280 1.608166 2.633244 3.621200 2.507290 3.988309
##  [8] 4.721993 3.776069 7.126221

t.test(x - y)

## 
##  One Sample t-test
## 
## data:  x - y
## t = -0.69081, df = 9, p-value = 0.5071
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.4789413  0.7869793
## sample estimates:
## mean of x 
## -0.345981

t.test(x, y, paired = TRUE)

## 
##  Paired t-test
## 
## data:  x and y
## t = -0.69081, df = 9, p-value = 0.5071
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.4789413  0.7869793
## sample estimates:
## mean of the differences 
##               -0.345981

x <- rlnorm(10, meanlog = 3, sdlog = 1)
y <- rlnorm(10, meanlog = 3.5, sdlog = 1)
list(x = x,y = y)

## $x
##  [1] 10.029806  7.141222 79.677275 23.761654  6.998416 16.235389 27.343574
##  [8] 24.912341  5.506661 12.352728
## 
## $y
##  [1]   9.149032 120.784095  29.724325  49.913021  77.843253  42.675081
##  [7]  38.263522  15.387308 260.455593 109.190266

wilcox.test(x - y)

## 
##  Wilcoxon signed rank test
## 
## data:  x - y
## V = 9, p-value = 0.06445
## alternative hypothesis: true location is not equal to 0

wilcox.test(x, y, paired = TRUE)

## 
##  Wilcoxon signed rank test
## 
## data:  x and y
## V = 9, p-value = 0.06445
## alternative hypothesis: true location shift is not equal to 0

Si los datos no están asociados, entonces

• si ambas muestras están normalmente distribuidas, o el conjunto de datos –escalados y centrados en cada muestra– está normalmente distribuido, se usa la prueba de t. En R, se usa la función t.test; la opción var.equal permite escoger entre la prueba calculada con la varianza de los datos, y la aproximación de Welch, que se usa cuando la varianzas no son iguales.
• si los datos no están normalmente distribuidos se usa la prueba U de Mann-Whitney –wilcox.test en R–.

x <- rnorm(10, mean = 3, sd = 1)
y <- rnorm(14, mean = 3.5, sd = 2)
list(x = x,y = y)

## $x
##  [1] 3.966946 4.299946 3.341676 2.564442 3.300222 2.965613 2.086873
##  [8] 3.419466 3.124921 2.672670
## 
## $y
##  [1] -2.0898613  5.9415309  2.1861442  5.9859275  3.1561266  0.7314367
##  [7]  1.6177214  1.7548402  1.8578332  0.4620650  4.1497187  0.7134045
## [13]  2.5924116  5.3051376

t.test(x, y, var.equal = FALSE)

## 
##  Welch Two Sample t-test
## 
## data:  x and y
## t = 1.1125, df = 15.833, p-value = 0.2825
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6529183  2.0922681
## sample estimates:
## mean of x mean of y 
##  3.174277  2.454603

x <- rlnorm(10, meanlog = 3, sdlog = 1)
y <- rlnorm(14, meanlog = 3.5, sdlog = 2)
list(x = x,y = y)

## $x
##  [1]  29.653909  23.152165 111.874130   9.743348  58.153463  52.115569
##  [7]  55.746236  30.161096  82.766524  30.055187
## 
## $y
##  [1]  216.565053    5.803752   12.543098  746.357529  108.242239
##  [6]    2.158818    0.664829    7.535584 3575.733941  154.265032
## [11]   21.202641    3.565533   26.496414  154.990273

wilcox.test(x, y)

## 
##  Wilcoxon rank sum test
## 
## data:  x and y
## W = 77, p-value = 0.7088
## alternative hypothesis: true location shift is not equal to 0

Para la comparación de tres o más grupos, las técnicas más habituales son

• si los datos de los distintos conjuntos de datos están normalmente distribuidos y presentan dispersiones comparables –homocedasticidad–, la técnica recomendada es el análisis de la varianza –aov en R–;
• si los datos no están normalmente distribuidos pero las distribuciones tienen formas parecidas y se mantiene la homocedaticidad, se usa la prueba de Kruskal-Wallis –en R, krukal.test–.

x <- c(rnorm(10, mean = 3, sd = 1),
       rnorm(8, mean = 3.2, sd = 1),
       rnorm(10, mean = 4, sd = 1))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))
test <- aov(x ~ g)
summary(test)

##             Df Sum Sq Mean Sq F value Pr(>F)  
## g            2  7.523   3.761   3.449 0.0476 *
## Residuals   25 27.268   1.091                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

x <- c(runif(10, min = 2, max = 4),
       runif(8, min = 1.5, max = 3.5),
       runif(10, min = 3, max = 5))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))
kruskal.test(x ~ g)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  x by g
## Kruskal-Wallis chi-squared = 14.279, df = 2, p-value = 0.0007932

Tanto el análisis de la varianza como la prueba de Kruskal-Wallis solo permiten contrastar si algún par de grupos presentan valores de localización distintos, pero no permiten identificar cuál o cuáles.

Para identificar qué pares de grupos presentan diferencias estadísticamente significativos se usan pruebas post-hoc

• En el caso del análisi de varianza, se puede usar la prueba de Tukey –TukeyHSD en R,
• Para la prueba de Kruskal-Wallis, se puede usar la prueba de Dunn –dunn.test en el paquete dunn.test de R–

Gráficamente y aunque con menos objetividad, la comparación de la loclización puede hacerse usando gráficos de caja.

x <- c(rnorm(10, mean = 3, sd = 1),
       rnorm(8, mean = 3.2, sd = 1),
       rnorm(10, mean = 4, sd = 1))
g <- factor(c(rep(1,10), rep(2,8), rep(3,10)))

xs <- data.frame(y = x, group = g)
  
ggplot(xs, aes(x=group,y=y)) + geom_boxplot() +
  coord_flip() + theme_bw()

Ajuste de modelos

El ajuste de modelos se expresa de forma habitual en R mediante un objeto formula.

En esta notación, se indica la relación entre variables mediante una expresión de tres términos donde la variable dependiente se indica a la izquierda y las variables dependientes a la derecha.

form1 <- y ~ x              # recta
class(form1)

## [1] "formula"

form1 <- y ~ log(x)         # logaritmo
form1 <- y ~ poly(x,4)      # polinomio de grado 4
form1 <- y ~ x + 0          # recta que pasa por el origen
form1 <- y ~ I(x^.5)        # raíz cuadrada

Ajuste de modelos – regresión lineal

El ajuste de un modelo lineal por mínimos cuadrados ordinarios se hace en R mediante la función lm.

fit1 <- lm(mpg ~ wt, data=mtcars)
summary(fit1)

## 
## Call:
## lm(formula = mpg ~ wt, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.5432 -2.3647 -0.1252  1.4096  6.8727 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
## wt           -5.3445     0.5591  -9.559 1.29e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.046 on 30 degrees of freedom
## Multiple R-squared:  0.7528, Adjusted R-squared:  0.7446 
## F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

Una vez ajustado el modelo, este puede usarse para predecir nuevos datos (o calcular los valores ajustados) mediante la función predict.

df <- data.frame(x=mtcars$wt, y=mtcars$mpg,
  predict(fit1, newdata=mtcars, interval="prediction"))
head(df)

##                       x    y      fit      lwr      upr
## Mazda RX4         2.620 21.0 23.28261 16.92894 29.63628
## Mazda RX4 Wag     2.875 21.0 21.91977 15.59072 28.24882
## Datsun 710        2.320 22.8 24.88595 18.48644 31.28546
## Hornet 4 Drive    3.215 21.4 20.10265 13.78568 26.41962
## Hornet Sportabout 3.440 18.7 18.90014 12.57806 25.22223
## Valiant           3.460 18.1 18.79325 12.47021 25.11630

El modelo se puede visualizar gráficamente a partir de estas predicciones.

df2 <- data.frame(
  wt = seq(min(mtcars$wt), max(mtcars$wt), length.out = 201),
  predict(fit1, 
  newdata = data.frame(wt = seq(min(mtcars$wt), max(mtcars$wt),
                             length.out = 201)), 
  interval="prediction"))

ggplot(mtcars, aes(x=wt, y=mpg)) +
  geom_ribbon(data=df2, mapping = aes(x=wt, ymin=lwr,ymax=upr),
              inherit.aes = FALSE, fill="#dddddd", alpha=0.8)+
  geom_point(shape=3) +
  geom_line(data=df2, mapping = aes(x=wt, y=fit)) +
  theme_bw()

Ajuste de modelos – comprobación

Para comprobar la adecuación del modelo y del método de ajuste usado (OLS) deben analizarse los residuales.

Estos deben ser aleatorios, estar normalmente distribuidos, mostrar homocedasticidad y no mostrar evidencia de puntos especialmente influyentes.

Estas comprobaciones pueden hacerse a partir de los datos de los residuales y las pruebas estadísticas oportunas.

fit1$residuals

##           Mazda RX4       Mazda RX4 Wag          Datsun 710 
##          -2.2826106          -0.9197704          -2.0859521 
##      Hornet 4 Drive   Hornet Sportabout             Valiant 
##           1.2973499          -0.2001440          -0.6932545 
##          Duster 360           Merc 240D            Merc 230 
##          -3.9053627           4.1637381           2.3499593 
##            Merc 280           Merc 280C          Merc 450SE 
##           0.2998560          -1.1001440           0.8668731 
##          Merc 450SL         Merc 450SLC  Cadillac Fleetwood 
##          -0.0502472          -1.8830236           1.1733496 
## Lincoln Continental   Chrysler Imperial            Fiat 128 
##           2.1032876           5.9810744           6.8727113 
##         Honda Civic      Toyota Corolla       Toyota Corona 
##           1.7461954           6.4219792          -2.6110037 
##    Dodge Challenger         AMC Javelin          Camaro Z28 
##          -2.9725862          -3.7268663          -3.4623553 
##    Pontiac Firebird           Fiat X1-9       Porsche 914-2 
##           2.4643670           0.3564263           0.1520430 
##        Lotus Europa      Ford Pantera L        Ferrari Dino 
##           1.2010593          -4.5431513          -2.7809399 
##       Maserati Bora          Volvo 142E 
##          -3.2053627          -1.0274952

Gráficamente, los mismas comprobaciones pueden hacerse a partir de la representación gráfica del modelo.

plot(fit1, which=1:6)

Estas mismas ideas son la base de la aplicación en R de las técnicas de ajuste multilineal (lm), ajuste lineal generalizado (glm) y ajuste no lineal por mínimos cuadrados (nls).